par(mfrow=c(1,4)) par(lend=1) ?dnbinom # Fix mu, vary theta plot(0:20, dnbinom(0:20,mu=1,size=0.1),type='h',lwd=4,xlab='count category',ylab='probability',ylim=c(0,.8)) plot(0:20, dnbinom(0:20,mu=1,size=0.1),type='h',lwd=4,xlab='count category',ylab='probability',ylim=c(0,.8)) mtext(side=3,line=.5,expression(list(mu==1,theta==0.1))) plot(0:20, dnbinom(0:20,mu=1,size=1),type='h',lwd=4,xlab='count category',ylab='probability',ylim=c(0,.8)) mtext(side=3,line=.5,expression(list(mu==1,theta==1))) plot(0:20, dnbinom(0:20,mu=1,size=10),type='h',lwd=4,xlab='count category',ylab='probability',ylim=c(0,.8)) mtext(side=3,line=.5,expression(list(mu==1,theta==10))) plot(0:20, dnbinom(0:20,mu=1,size=100),type='h',lwd=4,xlab='count category',ylab='probability',ylim=c(0,.8)) mtext(side=3,line=.5,expression(list(mu==1,theta==100))) # compare Poisson and NB par(mfrow=c(1,2)) plot(0:12, dnbinom(0:12,mu=1,size=100),type='h',lwd=4,xlab='count category',ylab='probability',ylim=c(0,.8)) mtext(side=3,line=.5,expression(paste('NB: ',list(mu==1,theta==100)))) plot(0:12, dpois(0:12,lambda=1), type='h',lwd=4,xlab='count category',ylab='probability',ylim=c(0,.8)) mtext(side=3,line=.5,expression(paste('Poisson: ',mu==1))) # fix theta, vary mu par(mfrow=c(1,3)) plot(0:40, dnbinom(0:40,mu=1,size=0.1),type='h',lwd=4,xlab='count category',ylab='probability',ylim=c(0,.8)) mtext(side=3,line=.5,expression(list(mu==1,theta==.1))) plot(0:40, dnbinom(0:40,mu=10,size=0.1),type='h',lwd=4,xlab='count category',ylab='probability',ylim=c(0,.8)) mtext(side=3,line=.5,expression(list(mu==10,theta==.1))) plot(0:40, dnbinom(0:40,mu=20,size=0.1),type='h',lwd=4,xlab='count category',ylab='probability',ylim=c(0,.8)) mtext(side=3,line=.5,expression(list(mu==20,theta==.1))) par(mfrow=c(1,1)) ##### Aphid data set ###### num.stems <- c(6,8,9,6,6,2,5,3,1,4) # data frame of tabulated data aphids <- data.frame(aphids=0:9, counts=num.stems) aphid.data <- rep(0:9,num.stems) poisson.LL <- function(lambda) sum(log(dpois(aphid.data, lambda))) poisson.negloglik <- function(lambda) -poisson.LL(lambda) out.pois <- nlm(poisson.negloglik,3) exp.pois <- c((dpois(0:8, out.pois$estimate)), 1-ppois(8,out.pois$estimate))*50 out.bar <- barplot(num.stems, ylim=c(0,11), names.arg=c(0:8,'9+')) points(out.bar, exp.pois, pch=16, cex=.9, type='o') legend('topright', 'Poisson model', pch=16, col=1, lty=1, cex=.9, bty='n') # version with two arguments NB.LL <- function(mu,theta) sum(log(dnbinom(aphid.data,mu=mu, size=theta))) NB.LL(2,3) # version with a vector argument NB.LL2 <- function(p) sum(log(dnbinom(aphid.data,mu=p[1], size=p[2]))) # need negative log-likelihood negNB.LL2 <- function(p) -sum(log(dnbinom(aphid.data,mu=p[1], size=p[2]))) out.nb <- nlm(negNB.LL2, c(2,3)) out.nb # predicted probabilities with tail category exp.p <- c(dnbinom(0:8,mu=out.nb$estimate[1], size=out.nb$estimate[2]), pnbinom(8,mu=out.nb$estimate[1], size=out.nb$estimate[2],lower.tail=F)) exp.p #expected counts exp.nb <- 50*exp.p #add NB predictions to graph points(out.bar, exp.nb, pch=15, cex=.9, type='o', col=2) exp.nb # collaps categories to pass the >5 rule exp.nb.short <- c(exp.nb[1:5], sum(exp.nb[6:7]), sum(exp.nb[8:10])) exp.nb.short # group observed values the same way obs.short <- c(num.stems[1:5], sum(num.stems[6:7]), sum(num.stems[8:10])) obs.short # Pearson statistic pearson <- sum((obs.short-exp.nb.short)^2/exp.nb.short) pearson # p-value 1-pchisq(pearson, length(obs.short)-1-2) # parametric bootstrap chisq.test(num.stems, p=exp.p, simulate.p.value=T, B=999) # fit negative binomial using glm.nb library(MASS) out.glmnb <- glm.nb(aphid.data~1) summary(out.glmnb) coef(out.glmnb) # undo log link exp(coef(out.glmnb)) names(out.glmnb) ### slugs data set #### slugs <- read.delim('ecol 563/slugsurvey.txt') slugs[1:8,] slug.table <- table(slugs$slugs, slugs$field) slug.table # bar plots of counts #stacked bar barplot(slug.table) # separate distributions barplot(slug.table, beside=T) coord1 <- barplot(slug.table, beside=T)->coord1 coord1 # add labels to bars coord1 <- barplot(slug.table, beside=T, names.arg=c(0:10,0:10), cex.names=0.75) # add labels for groups mtext(side=1, line=2.5, at=coord1[6,], text=c('Nursery','Rookery'), cex=.9) # interweaved bars barplot(t(slug.table), beside=T, legend=T) # lattice plot library(lattice) # make a data frame out of table slugtable <- data.frame(slug.table) # count categories are a factor variable slugtable$Var1 # panels of bar plots barchart(Freq~Var1|Var2, data=slugtable) barchart(Freq~Var1|Var2, data=slugtable, xlab='count category', col='grey70') xyplot(Freq~Var1|Var2, data=slugtable, xlab='count category', col='grey70', type='h', lwd=6) xyplot(Freq~Var1|Var2, data=slugtable, xlab='count category', col='grey70', type='h', lwd=8, lineend=1) # Single Poisson mean model negpoisLL<-function(p) { LL <- sum(log(dpois(slugs$slugs,lambda=p))) -LL} mean(slugs$slugs) negpoisLL(2) negpoisLL(1.775) out.pois1 <- nlm(negpoisLL,3) out.pois1 out.pois1 <- nlm(negpoisLL,1.8)->out.pois # separate mean Poisson model negpoisLL2<-function(p) { LL <- sum(log(dpois(slugs$slugs, lambda=p[1] + p[2]*(slugs$field=='Rookery')))) -LL} negpoisLL2(c(1,2)) out.pois2 <- nlm(negpoisLL2,c(1,2)) out.pois2 # likelihood ratio test that means are different 2*(out.pois1$minimum-out.pois2$minimum) LR <- 2*(out.pois1$minimum-out.pois2$minimum) 1-pchisq(LR,1) # fit same models using glm out1 <- glm(slugs~1, data=slugs ,family=poisson) out2 <- glm(slugs~field, data=slugs, family=poisson) # likelihood ratio test comparing models anova(out1, out2, test='Chisq') # sequential likelihood ratio tests anova(out2, test='Chisq') # Wald test summary(out2) # negative binomial version of the models out1.nb <- glm.nb(slugs~1, data=slugs) out2.nb <- glm.nb(slugs~field, data=slugs) anova(out1.nb, out2.nb) summary(out2.nb) # obtain predicted means predict(out2, type='response') # get means for each field type slugtable$mu <- predict(out2, type='response', newdata=data.frame(field=slugtable$Var2)) slugtable # calcualte probabilities of each category slugtable$p <- dpois(as.numeric(as.character(slugtable$Var1)), lambda=slugtable$mu) # add tail probability to the last category slugtable$p2 <- dpois(as.numeric(as.character(slugtable$Var1)), lambda=slugtable$mu)+ (slugtable$Var1==10)*ppois(as.numeric(as.character(slugtable$Var1)), lambda=slugtable$mu, lower.tail=F) sum(slugtable$p[1:11]) sum(slugtable$p2[1:11]) # obtain expected counts table(slugs$field) slugtable$exp <- 40*slugtable$p2 slugtable # add Poisson predictions to the graph xyplot(Freq~Var1|Var2, data=slugtable, xlab='count category', panel=function(x,y,subscripts) { panel.xyplot(x, y, type='h', col='grey70', lwd=8, lineend=1) panel.points(x, slugtable$exp[subscripts], pch=16, type='o')})